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Finite density simulations require dynamical fermions which are computationally demanding. 
We employ Ferrenberg and Swendsen reweighting method to reduce the number of ensembles 
needed. We use their method to do a combined reweighting in both j3 and quark number k on an 
ensemble generated by the Kentucky group. 
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1. Introduction 

Finite density simulations require dynamical fermions which are computationally demanding. 
When studying phase transitions we have to scan very finely both the temperature and density 
directions; a large number of ensembles need to be generated. Ferrenberg and Swendsen [[Tj] [||] 
showed that we can employ reweighting to reduce the number of ensembles needed. Using their 
methods we perform an extrapolation in j3 and then we do a combined reweighting in both j8 
and the quark number k. We take the ensembles generated at k = and k = 3, and used them to 
extrapolate to k = 6 case. The results of reweighting are also compared against the ones of that 
direct measurements. 



2. Simulation Background 

Finite density ensembles were generated using the canonical approach [J3j] - To build canonical 
partition function, we start from fugacity expansion of the grand canonical partition function, 

Z(V,T,li) = J £z c (V,T,k)e^ T , (2.1) 

k 

where k is the net quark number, Zq is the canonical partition function of the system. On a lattice, 
we can easily compute the Fourier transform of the grand canonical partition function 

Z(V,T,n) = J ' mJ@ijir9}ife- s *W-SfQiV,9,Y) (2.2) 
to get the canonical partition function 

1 f^ n 

Z c (V,T,k) = —J o &<!>e- ik *Z{V,T,n)\^ T . (2.3) 

We will specialize to the case of two degenerate flavors. After integrating out the fermion part, we 
get a simple expression 

Z c (V,T,k)= J ^Ve- s ^ v) det k M 2 (U), (2.4) 

where 

1 t 1% 

det k M 2 {U) = — d$e-' k UetM(m,n;U) 2 \v=i$T (2.5) 
27T Jo 

is the determinant projection to the fixed net quark number k. Then Zq is used to generate the 
ensembles at particular k. 



Simulations based on the action give in (2.4) cannot be carried out directly since the integrand 



is not positive definite. The determinant is split in two parts: a positive definite part that is used 
to generate the ensemble and a phase. For a more detailed presentation we refer the reader to the 
original paper [S]. 

The absolute value of the Polyakov loop is given by 

(\P\) = (2.6) 
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where 

l2lTl\ 1 N-l 



a(U) = , d^M 2 (£/) = I £ e-^detM(^) 2 (2.7) 



RedetjtM 2 ^) 



is the phase and () stands for the average over the ensemble generated with measure Redet£M 2 (£/) 
Quark chemical potential is defined by the expression below: 

N-l 

j 
1 

More relevant to our study is the baryon chemical potential: 



j N-l 
X — 



HB(n B ) = ~]n(e-*% nB . (2.9) 



P 



To test the reweighting method, a subset of the original ensembles was employed and the 
results of the extrapolation were compared against the results of the direct measurements. 

3. Algorithm 

The original reweighting method [jl]] employed an ensemble generated at a particular point 
in the parameter space to build an "induced" ensemble for a different value of parameters. The 
method is limited to a neighborhood of the original point in the parameter space. At distant points 
the extrapolation becomes unreliable due to poor overlap between the original and the target distri- 
bution. 

The multi-histogram method [0] addresses the limitations of the original method by employ- 
ing several ensembles. A carefully chosen set of ensembles, covering the area of interest in the 
parameter space, allows us to interpolate reliably. Thus, we can finely scan the area of interest with 
only a small set of ensembles. 

The basic idea behind the reweighting method is to use the histogram generated by the sim- 
ulated ensembles as an approximation of the probability distribution; this approximation is then 
used to approximate the probability distribution at a different point in the parameter space. Once 
we have an approximation for the probability distribution we can compute all observables. 

We used these methods in two scenarios: 

1. Fix k, vary j8. 

• Single histogram method: The histogram of only one ensemble is used in this case. 

• Multi-histogram method: The histograms of several ensembles generated at different 
points are combined. 

2. Vary both j3 and k. (Only Multi-histogram method is used.) 
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In our first case (fix k, vary j3), we approximate the "weight" of each configuration U by 

e s P (u) 

p ( u >P) = 7H s (in f 

This approximation allows us to compute averages at arbitrary values of j3. The input for this 
formula consists of R ensembles generated by weights given by Sp [ , . . . ,Sp R with n\ , . . . ,hr config- 
urations in each ensemble. 

To use the formula above we need to determine the parameters /„ — they are the "free ener- 
gies" of each ensemble n, 

exp{f n )=Y,P{UuPn) (3-2) 



where we sum over all the configurations Uj in the ensemble. Eq. (gjj) and ( p.2| ) are used to 
determine the free energy /„. We start with a guess for /„ and plug this into Eq. (3T) to get 
an approximation for the probability distribution P(U,p). Then, we take this approximation for 
P(U,P) and plug it into Eq. ( |3.2| ) to get a more refined guess of /„. We repeat this process until /„ 
converges. 

After convergence is reached, we can measure the observables at different /3 by 

For our second case (both j8 and k are varied), the "weight" for different configurations at 
different k and j3 should be calculated as 

P(U,k',P') = —^ „ .... . (3.4 

Lf= i m I Re det^.M 2 (U) k~ s ft {u) - fl 

exp(-fi) = y\P(Uj,ki,f5i) We sum over all configurations (3.5) 
The chemical potential can be measured by using the "weight" P(U,k,fi): 

j N-l 



X — 

N J=o 



£ e-^e-^AtXM 2 ^^) (3.6) 



Optimized convergence 



In our calculations, we found that the iteration presented above converges very slowly. For 
a more efficient convergence we use the following procedure: we know that (the k + X th 

iteration step) can be expressed as a function of /w (k th step) 

f(k+l) _ Q^(k) j Q is the yector function of f(k) (37) 



In our case, / is the free energy, G is the recurrent equation (3.5) for the free energy calculation. 
The fixed point has the property that /M = G(f^). At step k our approximation is /w, we 
write then: 

ft*) + 8f® = G{f {k) + 8f {k) ) where 8f [k) = f H - f® . (3.8) 
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Assuming that is close to then 5/W is small and we can approximate 

G(f [k) + 8f {k) ) » G(f {k) ) + VG • S/ w . (3.9) 

Using this we get: 

fl k) + 8tf^ Gi m + ^8ff (3.10) 
is the element of matrix = 

/f + S/f « GttfM) +M\fdff (3.11) 

so 

(/-M«) i7 5/f « G,-(/«) -/f (3.12) 

5/f ~ (Gitf®)-^) ■ (I-M^)-\j (3.13) 

We expect the matrix (/ — ) to be almost singular when we are approaching fixed point, 
so we use SVD to compute the matrix inversion. Our next step in the iteration is then p k+l ) = 

+ 5/W rather than G(/W). Using this method the convergence of the iteration is greatly 
accelerated. 



4. Results and Conclusion 

We determine the Polyakov loop using the single histogram method and the multi-histogram 
method and compare the results with the direct measurements. The results are plotted as a function 
of j6 , the gauge coupling constant. 
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Figure 1: Polyakov loop (Single histogram at fixed k = 3). All dots represent the ensembles where direct 
measurements were carried out. The ones surrounded with polygon are used as our reweighting input; the 
gray ones are not. The line represents the result of the extrapolation and the band represents the statistical 
errors. 

From single histogram reweighting of Polyakov loop, we can see from Fig. 1 the method is 
limited to a neighborhood of the original point j8 = 5.20 in the parameter space. The extrapolation 
at distant points becomes unreliable due to poor overlap. The multi-histogram method addresses 
this limitation, and by employing several ensembles, we expect to interpolate results reliably. 
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We start by employing three ensembles in Fig. 2. Two of them at /3 = 4.8 and j3 = 5.6 are the 
bounds of the region of interest and the third one is in the middle at j8 = 5.2. We plot the results 
in Fig. 2a; we see that the error bars are quite large between j8 = 4.9 and 5.1. When we add one 
more ensemble in the "problem" region, the plot changes significantly as we can see from Fig. 2b. 
The largest error bars are around j3 =5.3 — 5.4 now. Once we add a point at jS = 5.3 the error 
bars are reduced as we see from Fig. 2c. From this plot we also see that with only half of the 
generated ensembles the reweighting curve describes very well the data in the region of interest. 
From Fig. 2d, we see that adding more ensembles (total 10 ensembles) doesn't change the curve 
significantly — the only difference is that the error band gets smaller. 



Polyakov loop chemical potential 




Figure 3: We used ensembles generated at k = and k = 3 to extrapolate k = 6 case, and the plot is compared 
to the direct measurements. 
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We have also performed a combined reweighting in /3 and k in Fig. 3. We employed ten 
ensembles at k = and ten ensembles at k = 3 to extrapolate to k = 6. We used this method 
to compute both the Polyakov loop and the chemical potential. As we can see from Fig. 3 the 
extrapolated curve matches rather well the direct measurements. 

We find that multi-histogram method can be used effectively to reduce the number of ensem- 
bles needed for our simulations: for our j8 interpolations we find that we could use as little as half 
of the ensembles generated to get a much better scan of the temperature range. The way we envi- 
sion using this method is to generate few ensembles of first scanning the temperature range rather 
coarsely and then add some more ensembles at the points where "features" like sharp transitions 
appear in our plots. The process stops when adding new ensembles does not alter our plots signifi- 
cantly. We also find that interpolations in both temperature and quark number are feasible and we 
plan to employ these methods in our further studies. 
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